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ABSTRACT 

The large number of microlensing events discovered towards the Galactic Bulge bears 
the promise to reconstruct the stellar mass function. The more interesting issue con- 
cerning the mass function is certainly to probe its low mass end, near the region 
occupied by the brown dwarfs. However due to the source confusion, even if the distri- 
bution and the kinematics of the lenses are known, the estimation of the mass function 
is extremely biased at low masses. The blending due to the source confusion biases 
the duration of the event, which in turn dramatically biases the estimation of the 
mass of the lens. To overcome this problem we propose to use differential photometry 
of the microlensing events obtained using the image subtraction method. Differential 
photometry is free of any bias due to blending, however the drawback of differential 
photometry is that the baseline flux is unknown. In this paper we will show that even 
without knowing the baseline flux, purely differential photometry allow to estimate 
the mass function without any biases. The basis of the method is that taking the scalar 
product of the microlensing light curves with a given function and taking its sum over 
all the microlensing events is equivalent to project the mass function on another func- 
tion. This method demonstrates that there is a direct correspondancy between the 
space of the observations and the space of the mass function. Concerning the function 
to use in order to project the observations, we show that the principal components 
of the light curves are an optimal set. We also demonstrate that there is no addi- 
tional information about the distribution of the scalar products of the data beyond 
their sum (first order moment) . Higher order moments are only linear combination of 
the first order moment. Thus the sum of the projection on the principal components 
contains all the information, and translate in an equal number of projection of the 
mass function with functions associated with the principal components. To illustrate 
the method we simulate data sets consistent with the microlensing experiments. By 
using 1000 of these simulations, we show that for instance the exponent of the mass 
function can be reconstructed without any biases. 
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I INTRODUCTION. 

The success of the microlensing experiments has been im- 
pressive, several collaborations have reported the detection 
of microlensing amplification of stars, OGLE (Udalski et al. 
2000), MACHO (Alcock et al. 1998), EROS (Derue et al. 
1999)), DUO (Alard & Guibert 1997), VATT (Uglesich et 
al. 1999). In particular, the recent release by the OGLE 

II collaboration of more than 200 microlensing events is of 
great interest. One of the obvious promises of such data sets 
is the possibility to explore the mass function near the low 
mass end. Although, the relation between the mass function 
and the microlensing observations is not straightforward. In 
the classical scheme of analysis, one must first estimate the 
duration of the events by fitting the theoretical light curves, 



then compute the lensing rates, and finally relate these lens- 
ing rates to the mass function. The trouble is that fitting the 
theoretical light curve to the microlensing data to obtain the 
duration of the event is a highly degenerated process (Alard 
1997, Han 1997, Wozniak & Paczynski 1997). Due to the 
source confusion in crowded fields, it is almost impossible to 
estimate the flux of the amplified source reliably. The flux 
of the source is severely biased by blending with neighbor- 
ing sources. In such case, an over-estimation of the baseline 
flux will result in an under-estimation of the duration of the 
event. This bias is severe for unresolved stars and will re- 
sult in a reduction in the estimation of the duration of a 
factor of 2,3, or even more. Since the bias on the mass of 
the lens goes like the square of the bias on the duration, the 
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relevant bias on the mass function will be very large. It was 
already shown (Alard 1997, Han 1997) that the contribu- 
tion of these unresolved, highly blended stars to the lensing 
rated was dominant at short durations. The consequence is 
that the estimation of the mass function at low mass end 
will be largely over-estimated, suggesting the existence of a 
large number of brown dwarf that actually do not exists. 
It is obvious that as long as the bias due to the blended 
sources has not been solved, any correct estimation of the 
mass function is impossible. This paper proposes a solution 
to this problem. We will see that using differential photom- 
etry obtained using the image subtraction method (Alard 
& Lupton 1998, Alard 2000) it is possible to estimate the 
mass function, even without knowing the baseline flux of 
the source. This method assumes that the distribution and 
kinematics of the lenses are known with good accuracy. It 
is important to emphasize that the uncertainties related to 
the structure of the Galaxy in the bulge region are several 
order of magnitude smaller than the uncertainties due to the 
blending bias. Many good model of the kinematic and struc- 
ture of the bulge of our Galaxy are already available (Zhao, 
1996, Fux 1997, Bissantz et ai, 1997, Binney, Gerhard, & 
Spergel 1997) 



2 THE METHOD. 
2.1 Introduction 

When analyzing microlensing observations one has to deal 
with the light curves of a number N of microlensing can- 
didates. We will assume that we have purely differential 
photometry only, with a baseline flux equal to zero. We em- 
phasize that high quality differential photometry can always 
be obtained after the events have been detected by classi- 
cal methods by using the image subtraction method (Alard 
1999). These light curves will be represented by the symbol 
Si(t), (i = 1, -, N). Assuming that un-amplified flux of the 
source is A, the expression of the Si(t) as a function of time 
will be given by the theoretical microlensing amplification 
formula: 



Si{t) 



With: 



{A- 



u{tf + 2 



(i) 



«(*) 



<+ — 
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(we recall that we are working in the assumption that we 
have at our disposal a differential flux only, with for conve- 
nience a baseline flux equal to 0). 

The most important issue is of course how do we relate these 
observations to the microlensing parameters: the impact pa- 
rameter, uo, and the time to cross the Einstein ring: tE. And 
furthermore how do we relate the microlensing parameters 
to the mass function itself. The first serious difficulty we 
encounter is that of course uo,tE, cannot be extracted eas- 
ily by fitting the light curve due to the extreme blending 
degeneracy (Alard 1997, Han 1997, Wozniak & Paczynski 
1997). Thus our first and most essential step will be to de- 
rive a better set of parameters. This set of parameters will 



have to be non degenerated with respect to blending. The 
most simple and most natural way to deal with such prob- 
lem is to express the observations as a linear combination 
of a small number of vectors. The best way to make such 
a decomposition is known as principal components analysis, 
it is equivalent to an eigen value decomposition. Thus all 
we have to do is to look for the firsts eigen functions, and 
to calculate the scalar product of the components with our 
vectors (time series) of observations. 



2.2 Normalization 

Before making the principal components analysis we have 
to consider that one of the parameters, the amplitude of the 
magnified sources does not contain any information about 
the lensing event. Actually the reason of the the degener- 
acy of the fit is precisely the unknown amplitude of the 
source. The only relevant and meaningful parameters are 
uo and tE - Thus it is important to derive an estimator which 
is independent of the amplitude. Since the amplitude is a 
linear parameter, it is sufficient to normalize the vector of 
observations in order to get rid of the amplitude. This nor- 
malization can be made in the sense of the scalar product 
we are going to use to make our principal components de- 
composition. In such case, the normalized light curve Si(t) 
can be expressed as: 



Si(t) = - x 



u{tf + 2 



a U(t) y/u(t) 2 +4 



And: 
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Note that A(t) does not depend any more on the amplitude, 
but only on u and ts- 

The most simple way to apply the previous normalization 
to the observations Sj is to calculate directly the modulus 
of \Sj\ from the observations. Although it is obvious that 
calculating directly \Si\ from the observations may not be 
optimal. To estimate \S%\ we will use the following approach: 
first we fit a microlensing model to the observations. Even 
if the fit is very degenerated, one can always find a good so- 
lution which fits the data well. We will call this fitted model 
Si. Once we get this light curve solution we estimate using 
the following formula: 



15, 



Si(t] dt 



We will use the same approach for all our other calculations, 
in the same way we will estimate the principal components 
not by taking directly the data, but by taking the models we 
have fitted to the light curves. And finally the scalar prod- 
ucts with the principal components will be also estimated 
by cross products with the theoretical models. 
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2.3 Principal components 



following formulae: 



As we have already explain, we will first fit a microlensing 
model Si to the light curves, normalize this series of N vec- 
tors to obtain the series of vectors Si, and then calculate the 
principal components of these Si vectors. Usually it is pos- 
sible to express almost all the information with a reduced 
number of principal components with the following linear 
decomposition: 

S i = Y J a 'V Pi® 

j 

Where ay is the projection of time series number i, 
Si (uo,tE,t) on the principal component number j: 

Oij (u ,t E ) = J Si (u ,t E ,t) Pj (t) dt 

Note that ay is a function of uo and tE: ay = /3 (uo,ts) 

For instance a typical set of 100 microlensing light curves 
can be expressed as the combination of only 4 components 
with an accuracy good to 1 %. The other principal compo- 
nents contains very little additional information, but mostly 
noise. 



2.4 Statistical distribution of the projections on 
the principal components. 

We consider that most of the information is contained in 
a number Nc of principal components. Thus almost all the 
information can be extracted from the distribution of the 
projection of the Ni time series on the Nc principal compo- 
nents. The first meaningful quantity concerning the distribu- 
tion of the ay is the first order moment of the distribution: 



i ) — ^ ] a ij 



The second interesting quantity is the second order moment 
of the distribution: 



( a2 j) 



But it is important to notice, that ay = [f3 (uo,tE)] , can be 
written as a scalar product of the data with a given function. 
Since we can always find a function -F(t) such that: 



ay = [13 (u ,t E )] 



St (u ,t E ,t) F(t) dt 



(2) 



Proof: 

To prove that F(t) exists we have to show that Eq. (2) can 
be satisfayed for any point in the space of the parameters 
(uo, ts). We can map this parameter space by using a regular 
grid with a step as small as we like. This grid will have an 
almost infinite number of rows (N u ) and columns (Nt). The 
total number of points in the grid is: Ng = N u x Nt. We 
define F(t) by sampling this function in Ng points within 
the limits of the parameter space of the variable t. Since Ng 
is very large, one can always approximate Eq. (2) with the 
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The above equation can be written Ng times (for all the 
points in the grid). Since we have also Ng unknown val- 
ues F m , we can write a full system of Ng linear equations, 
which will allow to find the Ng values F m which represents 
the function F(t). Thus the function F(t) exists for all val- 
ues of uo and tE. 

Using the principal components decomposition of Si we can 
write: 



N C 

£ 



Pj (t) F(t) dt 



Thus, finally we see that the distribution of the second mo- 
ment is not more than a linear combination of the the Nc 
firsts moments (mean). Thus the mean of the second mo- 
ment will not bring any additional information with respect 
to the mean. With a similar reasoning it is possible to show 
that the same property is true for the N th moment. 

Consequently we can conclude that all the information con- 
cerning the distribution of the ay is contained in the mean 
of these components. No additional information (uncorre- 
cted information) will be found by looking at higher order 
moments. 



2.5 From the principal components to the mass 
function. 

It is possible to re-express the previous definition of (aj) 
by using the number density p(uo,tE, M) to observe a given 
amplification of parameters (uo,tE) with a lens of mass M. 
In such case, the sum can be approximated very closely with 
an integral expression: 



( a i) 



p(u ,t E ,M) aj (uo.ts) duo dt E dM 



If we define the efficiency function of the experiment, 
0(uo,tE), the lensing rates for 1 solar mass lenses, To(i) 
and the mass function, 4>(M) the formulae for p reads: 



p(u ,t E ,M) = To 



<t>(M) e(u ,t E ) 



Leading to the following expression for (aj) 

tE 



(aj) 



'M 



<t>{M) 9 (uo, t E ) a t (u , t E ) du dt B dM 



It is easy to re-arrange this integral in the following way: 
(o 3 -> ~ J 4>(M)<t>(M) dM 

With: 



ip(M) = 



To 



6(u ,t E ) ay (u ,t E ) du dt E (3) 



Thus we see that basically projecting the microlensing light 
curve on a basis of function and taking the statistical sum is 
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equivalent to projecting the mass function itself on another 
function. Consequently, we see that a projection in the space 
of the observation (the light curve) is directly equivalent to 
a projection in the space of the mass function. The problem 
of finding the optimal set of projections for the observations, 
is solved by using the principal component method. Then all 
we have to do is to calculate the "image" of these principal 
components in the space of the mass function. 



3 MONTE-CARLO SIMULATIONS. 
3.1 Introduction. 

To illustrate this method and show how the mass function 
can be reconstructed, we will use a series of Monte-Carlo 
simulations. We will simulate microlensing events by select- 
ing the parameters of the events according to the density 
distribution, p(uo,tB, M). One additional parameter we will 
have to select is the amplitude of the source. To get random 
variables which reproduces the distribution p we need to de- 
compose the problem. Basically p(uo,tE,M) has 3 parts: 



The rates, Vm^e) 



t E 



iM(t B ) = r (— ) 



For To we will adopt the analytical expression given in Mao 
& Paczyhski (1996). 

The efficiency function, 0(uo,tE): 

For the efficiency we will adopt the following criteria: an 
event is detected if it has a minimum number (N m ) of data 
points above a 5 a threshold. For a given duration tE this 
criteria is simply equivalent to put a threshold on the im- 
pact parameter uo. All we have to do is to search for the 
maximum value ut of uo such that the amplification of at 
least N m data points is above 5 a. If the sampling is even, 
with a time step St, it just mean that the amplification must 
be larger than 5 a in a window of time around the maxima 
AT = N m St. Provided the maxima is at the origin of the 
time axis, it will finally result in the condition: 



Sj I ut, tE 



AT\ 



> 5 a 



(4) 



The distribution of no itself is uniform, to get et set of 
wo values we will use a random generator which will pro- 
vide a uniform distribution between and the maximum 
value for uo, um($e)- um^e) corresponds to the impact pa- 
rameter for the brightest source given by Eq xx. Once we 
have selected tE from the distribution Tm^e), we will cal- 
culate uo^thresh using Eq xx, if our random uo value is above 
ua.thresh, (uo, ts) = 0, otherwise 6(uo,tE) = 1. To sum- 
marize: 



O (u ,t E ) = 
The mass function, (j>(M) 




In all our simulation we will adopt a pure power law ex- 
pression for the mass function, with a lower (Mi) and an 
upper cut-off (Ms). 

f M~ a Mi < M < M s 
4>(M) = <^ 

I otherwise 

The amplitude: 

we will assume that the the flux of the un-amplified source 
amplitude distribution is a power law with an exponent of -2 
(Zhao, et al. 1995), in a given range of amplitude A max and 
Amin. In our simulation, the amplitudes will be generated 
by Monte-Carlo method, using this power law. 

3.2 Description of the implementation. 

To implement our Monte-Carlo simulation we will select the 
random variables in the following order: 

- The amplitude of the source, according to the probabil- 
ity law A~ v 

- The mass of the lens, according to the probability law A~ a 

- The duration of the event tE using the distribution, Tm^e) 

- The impact parameter uo using an uniform distribution in 
the range, 0, mm (ie). 

- Once t e and uo are selected we apply the efficiency cut-off, 
using the function O (uo,tE). 

This procedure produces a complete set of variables, A, tE, 
uo for the events passing the efficiency cut-off. Form the 
variable we compute the microlensing light curve using Eq. 
1. We add noise to the light curves according to the Poisson 
statistics. Using this procedure, it is possible to simulate a 
set of microlensing light curves Sj. 



3.3 Example 

In this example we will take an amplitude range which is 
typical for microlensing images towards the Galactic Bulge: 



-4,; 



10 6 , A min = 10 



For the mass function we take the following parameters: 

f M~ 2 0.05 < M < 10.0 
</>(M) = I 

K otherwise 
Using these parameters we simulated several series of mi- 
crolensing light curves. For illustration, we extracted a sam- 
ple of light curves from one of these simulations, they are 
presented in Fig. 1. The parameters were adjusted in order 
that in a simulation the total number of microlensing events 
simulated be close to 100. Then we proceeded to the princi- 
pal components decomposition. The orthonormal set of prin- 
cipal components was calculated using a singular value de- 
composition. The first 4 principal components are presented 
in Fig. 2. 

The functions ip(M) corresponding to the principal 
components in the space of the mass function can be com- 
puted using Eq. 3. Projecting the data on the princi- 
pal components is equivalent to projecting the mass func- 
tion on ip(M). To illustrate our discussion, an example of 
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Figure 1. Example of light curves as obtained using the Monte-Carlo simulation procedure described above. 



ip(M) functions calculated using the settings of the previ- 
ous Monte-Carlo simulation are presented in Fig. 3. 



3.4 Estimating the mass function. 

To illustrate the ability of the method to reconstruct the 
mass function we will assume that exponent of the mass the 
mass function, a is unknown. For each simulation we will 
calculate the 4 principal component of the set of microlens- 
ing light curves we have simulated. We will compute the sum 
of the projection on a principal component of all the light 



curves < aj >. Using Eq. 3 we will calculate the function 
which corresponds to each principal components in the space 
of the mass function. Then all we have to do is to compare 
the projection of a trial mass function with < aj >. The trial 
mass function which match the 4 < aj > as close as possible 
(in a least-square tens) will be the best mass function. This 
procedure was applied for 1000 simulation, each time the 
program derived the best value of the exponent of a. The 
histogram of the values of alpha for the 1000 simulations is 
presented in Fig. 4. 
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Figure 2. The first 4 principal components corresponding to the simulation presented in Fig. 1. 
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Figure 3. The function tp(M) which corresponds to the principal components. In the right panel we present the 4 original function, 
while in the other panel we present an orthonormal linear combination of these 4 functions. (Note that in the right panel the function 
were scaled to have the same amplitude). 
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Figure 4. Histogram of the values of Alpha (the exponent of the mass function) found by fitting the Monte-Carlo simulation by 
projecting the data in the space of the mass function. 
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